evolution_of_species¶Estas celdas es para importar la librería desde GitHub y que se mantenga actualizado. Hay que descomentar las líneas si se utiliza desde google colab.
# !git clone https://github.com/asolnn2a8/evolution_of_species.git
# %cd evolution_of_species
# !git pull origin master
# !pip install -r requirements.txt
# !python setup.py install
# %cd ..
import numpy as np
from matplotlib import pyplot as plt
from perthame_pde.model.solver import *
from perthame_pde.plotters import *
# Función para calcular cuanto se demora una función
import time
def timeit(method):
def timed(*args, **kw):
ts = time.time()
result = method(*args, **kw)
te = time.time()
print('%2.2f ms' % ((te - ts) * 1000))
return result
return timed
Se definen las funciones \begin{equation} F(t, x) = \int K(x, y) R(t, y) dy \end{equation} y \begin{equation} G(t, y) = \int r(x) K(x, y) u(t, x) dx \end{equation}
El proposito de las clases IntegrarF e IntegrarG es proporcionar una forma sencilla de realizar estos cálculos de forma eficiente.
Para esto representa estas funciones por medio de matrices, asumiendo que $K$, $R$, $r$ y $u$ vienen en forma matricial. Además de que una vez fijadas $K$ o $r$, luego no será necesario llamarlas en cada instante; sólo será necesario actualizar $R$ y $u$.
La forma de utilizar es definiendo instancias de la forma
F = FunctionarF(K, R_0, x_lims=(0, 1), y_lims=(0, 1), dt=0.01, N=100, M=100, T=100)
G = FunctionarG(r, K, u_0, **PARAMETROS_DISCRETIZACION)
donde los parámetros de discretización son opcionales
Para calcular los valores basta con preguntar la componente $(n, j)$-ésima y lo cálcula al mismo tiempo que lo entrega. También se pueden actualizar las matrices $R$ y $u$ en una misma línea. Además, se puede consultar por la matriz que lo representa. Por ejemplo:
# Retorna el valor de una componente
F[3, 4]
# Retorna el valor de una fila
G[3]
# Actualiza la matriz R y retorna la fila
F.actualize_row(R[3], 3)[3]
# Consulta por su matriz
G.matrix
from perthame_pde.model.integral_functional import FunctionalF, FunctionalG
def f1(x, y): return y ** 2
def f2(t, x): return np.ones_like(x)
N, M, T = 100, 110, 120
x_lims, y_lims, dt = (0, 1), (0, 2), 0.01
x = np.linspace(*x_lims, N + 2)
y = np.linspace(*y_lims, M + 2)
t = np.array([dt * n for n in range(T + 1)])
BOUNDARIES_CONFIG = {
"x_lims": x_lims,
"y_lims": y_lims,
"dt": dt,
"N": N,
"M": M,
"T": T,
}
r = np.linspace(-1, 0, N + 2)
F1 = np.array([f1(x_, y) for x_ in x])
F2 = np.array([f2(t_, x) for t_ in t])
# print(transform_to_discrete_function(F1, "K", N+2, x, M+2, y))
G = FunctionalG(r, F1, F2[0], **BOUNDARIES_CONFIG)
G.actualize_row(F2[3], 3)
print(G[3, 4])
print(G.x.lims)
-0.0025971917863809755 (0, 1)
Dado que tenemos dos métodos para calcular las Ecs. de Perthame, en donde la única diferencia radica en la actualización de $R$, se diseñaran dos solver que encapsulen las distintas formas de calcularlo.
El primer método consiste en resolve la EDO \begin{equation} \partial_t R(t, y) = -m_2(y) R(t, y) + R_{in}(y) - R(t, y) G(t, y) \end{equation} con $G(t, y)$ de la forma que se comentó anteriormente.
El segundo método consiste en asignar a $R$ el siguiente valor: \begin{equation} R(t, y) = \frac{R_{in}(y)}{m_2(y) + G(t, y)} \end{equation}
El primer solver utiliza el siguiente esquema de actualización \begin{equation} R^{n+1} = (1 - \Delta t \cdot (m_2 + G^{n})) \odot R^{n} + \Delta t R_{in} \end{equation}
donde $\odot$ representa la multiplicación componente a componente.
El segundo solver utiliza el siguiente esquema de actualización \begin{equation} R^{n+1} = \frac{ R_{in} }{ m_2 + G^{n+1} } \end{equation}
Y para el tercer solver se utilizó una aproximación centrada de segundo orden para la derivada temporal, resultando en el siguiente esquema: \begin{equation} R^{n+1} = R^{n-1} - 2 \Delta t \cdot (m_2 + G^{n}) \odot R^{n} + 2 \Delta t R_{in} \end{equation}
Para ocuparlo basta con generar una instancia de alguno de los solvers:
solver = Solver1R(r, K, u_0, R_0)
donde los parámetros de actualización son opcionales.
Al igual que antes, se puede actualizar u. Además, se puede buscar la componente $(n, k)$-ésima de la matriz.
# Calcula y retorna la fila n
solver[n]
# Calcula y retorna la componente n, k
solver[n, k]
# Actualiza u y retorna la componente
solver.actualize_row(u[2], 2)[n, k]
# Actualiza el solver al paso n + 1
solver.actualize_step_np1(1)
# Calcula la masa total (en este caso, los recursos totales) a lo largo del tiempo
solver.calculate_total_mass()
# Retorna una función R(t, y) interpolada
R = solver.function_interpolated()
from perthame_pde.model.solver import Solver1R, Solver2R, Solver3R
La ecuación
$$ \partial_t u(t,x) = u(t,x)(-m_1(x) + r(x)F(t,x)) + \epsilon \Delta u(t,x) $$se puede reescribir usando $A(t,x) = -m_1(x) + r(x)F(t,x)$ y se ve de la siguiente forma:
$$ \partial_t u(t,x) - \epsilon \Delta u(t,x) = A(t,x) u(t,x) $$Que es muy similar a la ecuación del calor, salvo que en este caso $f(t,x) = A(t,x)u(t,x)$ depende de $u$, y el Laplaciano está siendo perturbado por un término $\epsilon$. Es por ello, que se adaptará el esquema de discretización de la ecuación del calor para resolver este problema, incorporando la dependencia de $u$ de la función $f$ y la pertubación.
Entonces, reemplazando en la ecuación con las diferencias finitas forward de primer y centrada de segundo orden, se obtiene:
$$ \frac{u^{n+1}_j - u^{n}_j}{\Delta t} - \epsilon \frac{u^{n}_{j-1} - 2u^{n}_j - u^{n}_{j-1}}{h^2} = A^n_j u^n_j $$que al despejar $u_j^{n+1}$ es de la siguiente forma:
$$ u_j^{n+1} = \alpha \ u_{j-1}^{n} + \beta_{j}^{n} \ u_{j}^{n} + \gamma \ u_{j+1}^{n} $$con $\alpha$, $\beta$ y $\gamma$ $$ \alpha = \epsilon \frac{\Delta t}{h^2} \qquad \beta_{j}^{n} = 1 - 2\alpha + \Delta t A^n_j \qquad \gamma = \alpha $$
from perthame_pde.model.solver import Solver1U, Solver2U
from perthame_pde.model.solver import solve_perthame
def norm_func(x, mu=0., sig=1.) -> float:
num = np.exp(-(x - mu)**2 / (2 * sig**2))
den = (sig * np.sqrt(2 * np.pi))
return num / den
def resolver_perthame_gaussiano(
m_1, m_2,
eps, sig_K, r,
M_in, mu_in, sig_in,
u_max, mu_u, sig_u,
R0_max, mu_R0, sig_R0,
solver_u=None, solver_R=None,
x_lims=(0, 1), y_lims=(0, 1), N=100, M=100, dt=0.01, T=100,
verbose=True, reports_every=10,
):
if solver_u is None:
solver_u = Solver2U
if solver_R is None:
solver_R = Solver1R
def _K(x, y, sig_K=sig_K) -> float:
return norm_func(x, mu=y, sig=sig_K)
def _R_in(y, M_in=M_in, mu_in=mu_in, sig_in=sig_in) -> float:
return M_in * norm_func(y, mu=mu_in, sig=sig_in)
def _m_1(x) -> float:
return m_1
def _m_2(y) -> float:
return m_2
def _r(x) -> float:
return r
# Distribuciones Gaussianas
def _u_0(x, u_max=u_max, mu_u=mu_u, sig_u=sig_u) -> float:
return u_max * norm_func(x, mu=mu_u, sig=sig_u)
def _R_0(y, R0_max=R0_max, mu_R0=mu_R0, sig_R0=sig_R0) -> float:
return R0_max * norm_func(y, mu=mu_R0, sig=sig_R0)
u, R = solve_perthame(
_u_0, _R_0, _r, _R_in, _m_1, _m_2, _K, eps,
solver_R=solver_R, solver_u=solver_u,
verbose=verbose, reports_every=reports_every,
x_lims=x_lims, y_lims=y_lims, N=N, M=M, dt=dt, T=T
)
return u, R
Con las hipótesis necesarias, si $\int \frac{R_{in}(y)}{m_2(y)}|\ln R^0(y)|dy < \infty$ y se tiene (para el caso Gaussiano) \begin{equation} \bar{m}_1 \bar{m}_2 \geq \frac{ r M_{in} }{ \sqrt{2\pi(\sigma_K^2 + \sigma_{in}^2)} } \end{equation} Entonces la solución se extinguirá, es decir, que $\int u(t, x) dx$ converge a cero mientras que $R(t, y)$ converge.
Se utilizaron $\varepsilon=0.001$, $\sigma_K=0.6$, $\sigma_{in}=1$, $r=1$, $M_{in}=3$ y $\bar{m}_1=\bar{m}_2=1.1$
u, R = resolver_perthame_gaussiano(
m_1=1.1, m_2=1.1,
eps=0.001, sig_K=0.5, r=1,
M_in=3, mu_in=0, sig_in=1,
u_max=0.6, mu_u=-0.3, sig_u=1,
R0_max=0.6, mu_R0=0, sig_R0=1,
solver_u=Solver2U, solver_R=Solver1R,
x_lims=(-4, 4), y_lims=(-4, 4), N=400, M=100, dt=0.01, T=1000
)
plot_total_mass(u, R)
plot_phase_plane(u, R, fig_config={"figsize": (12.5, 5)})
plot_contour_line(u, n_lines=300)
plt.show()
plot_anim(u, R, speed=5.0, t_min=None, t_max=None, ax_config={"ylim": (0, 1)})
Creating arrays... Building solvers... Starting iterations Iteration n = 999, Total Time: 21.036[seg], Time/Iter: 15.994[ms], avg. Time/Iter: 21.036[ms]
u, R = resolver_perthame_gaussiano(
m_1=1.1, m_2=1.1,
eps=0.001, sig_K=0.5, r=1,
M_in=3, mu_in=0, sig_in=1,
u_max=0.6, mu_u=-0.3, sig_u=1,
R0_max=0.6, mu_R0=0, sig_R0=1,
solver_u=Solver2U, solver_R=Solver2R,
x_lims=(-4, 4), y_lims=(-4, 4), N=400, M=100, dt=0.01, T=1000
)
plot_total_mass(u, R)
plot_phase_plane(u, R, fig_config={"figsize": (12.5, 5)})
plot_contour_line(u, n_lines=300)
plt.show()
plot_anim(u, R, speed=5.0, t_min=None, t_max=None, ax_config={"ylim": (0, 1)})
Creating arrays... Building solvers... Starting iterations Iteration n = 999, Total Time: 19.960[seg], Time/Iter: 16.000[ms], avg. Time/Iter: 19.960[ms]
Se utilizaron $\varepsilon=0.001$, $\sigma_K=0.6$, $\sigma_{in}=1$, $r=1$, $M_{in}=3$ y $\bar{m}_1=\bar{m}_2=1$
u, R = resolver_perthame_gaussiano(
m_1=1., m_2=1.,
eps=0.001, sig_K=0.5, r=1,
M_in=3, mu_in=0, sig_in=1,
u_max=0.6, mu_u=-0.3, sig_u=1,
R0_max=0.6, mu_R0=0, sig_R0=1,
solver_u=Solver2U, solver_R=Solver1R,
x_lims=(-4, 4), y_lims=(-4, 4), N=400, M=100, dt=0.01, T=1000
)
plot_total_mass(u, R)
plot_phase_plane(u, R, fig_config={"figsize": (12.5, 5)})
plot_contour_line(u, n_lines=300)
plt.show()
plot_anim(u, R, speed=5.0, t_min=None, t_max=None, ax_config={"ylim": (0, 1)})
Creating arrays... Building solvers... Starting iterations Iteration n = 999, Total Time: 17.017[seg], Time/Iter: 14.023[ms], avg. Time/Iter: 17.017[ms]
u, R = resolver_perthame_gaussiano(
m_1=1., m_2=1.,
eps=0.001, sig_K=0.5, r=1,
M_in=3, mu_in=0, sig_in=1,
u_max=0.6, mu_u=-0.3, sig_u=1,
R0_max=0.6, mu_R0=0, sig_R0=1,
solver_u=Solver2U, solver_R=Solver2R,
x_lims=(-4, 4), y_lims=(-4, 4), N=400, M=100, dt=0.01, T=1000
)
plot_total_mass(u, R)
plot_phase_plane(u, R, fig_config={"figsize": (12.5, 5)})
plot_contour_line(u, n_lines=300)
plt.show()
plot_anim(u, R, speed=5.0, t_min=None, t_max=None, ax_config={"ylim": (0, 1)})
Creating arrays... Building solvers... Starting iterations Iteration n = 999, Total Time: 17.352[seg], Time/Iter: 16.014[ms], avg. Time/Iter: 17.352[ms]
Se eligieron:
u, R = resolver_perthame_gaussiano(
m_1=1.5, m_2=0.5,
eps=0.001, sig_K=0.05, r=8,
M_in=5, mu_in=0, sig_in=2.5,
u_max=3, mu_u=-0.3, sig_u=0.01,
R0_max=5, mu_R0=0, sig_R0=1,
solver_u=Solver2U, solver_R=Solver1R,
x_lims=(-5, 5), y_lims=(-5, 5), N=500, M=150, dt=0.01, T=1000
)
plot_total_mass(u, R)
plot_phase_plane(u, R, fig_config={"figsize": (12.5, 5)})
plot_contour_line(u, n_lines=300)
plt.show()
plot_anim(u, R, speed=5.0, t_min=None, t_max=None, ax_config={"ylim": (0, 4)})
Creating arrays... Building solvers... Starting iterations Iteration n = 999, Total Time: 26.653[seg], Time/Iter: 35.995[ms], avg. Time/Iter: 26.653[ms]
u, R = resolver_perthame_gaussiano(
m_1=1.5, m_2=0.5,
eps=0.001, sig_K=0.05, r=8,
M_in=5, mu_in=0, sig_in=2.5,
u_max=3, mu_u=-0.3, sig_u=0.01,
R0_max=5, mu_R0=0, sig_R0=1,
solver_u=Solver2U, solver_R=Solver2R,
x_lims=(-5, 5), y_lims=(-5, 5), N=500, M=150, dt=0.01, T=1000
)
plot_total_mass(u, R)
plot_phase_plane(u, R, fig_config={"figsize": (12.5, 5)})
plot_contour_line(u, n_lines=300)
plt.show()
plot_anim(u, R, speed=5.0, t_min=None, t_max=None, ax_config={"ylim": (0, 4)})
Creating arrays... Building solvers... Starting iterations Iteration n = 999, Total Time: 25.038[seg], Time/Iter: 24.996[ms], avg. Time/Iter: 25.038[ms]
u, R = resolver_perthame_gaussiano(
m_1=0.5, m_2=0.5,
eps=0.001, sig_K=0.5, r=1,
M_in=3, mu_in=0, sig_in=3,
u_max=3.2, mu_u=-0.3, sig_u=0.1,
R0_max=0.6, mu_R0=0, sig_R0=1,
solver_u=Solver2U, solver_R=Solver1R,
x_lims=(-4, 4), y_lims=(-4, 4), N=400, M=100, dt=0.01, T=1000
)
plot_total_mass(u, R)
plot_phase_plane(u, R, fig_config={"figsize": (12.5, 5)})
plot_contour_line(u, n_lines=300)
plt.show()
plot_anim(u, R, speed=5.0, t_min=None, t_max=None, ax_config={"ylim": (0, 1)})
Creating arrays... Building solvers... Starting iterations Iteration n = 999, Total Time: 17.553[seg], Time/Iter: 20.975[ms], avg. Time/Iter: 17.553[ms]
u, R = resolver_perthame_gaussiano(
m_1=0.5, m_2=0.5,
eps=0.001, sig_K=0.5, r=1,
M_in=3, mu_in=0, sig_in=1,
u_max=3.2, mu_u=-0.3, sig_u=0.01,
R0_max=0.6, mu_R0=0, sig_R0=1,
solver_u=Solver2U, solver_R=Solver1R,
x_lims=(-4, 4), y_lims=(-4, 4), N=400, M=100, dt=0.01, T=1000
)
plot_total_mass(u, R)
plot_phase_plane(u, R, fig_config={"figsize": (12.5, 5)})
plot_contour_line(u, n_lines=300)
plt.show()
plot_anim(u, R, speed=5.0, t_min=None, t_max=None, ax_config={"ylim": (0, 4)})
Creating arrays... Building solvers... Starting iterations Iteration n = 999, Total Time: 19.230[seg], Time/Iter: 13.992[ms], avg. Time/Iter: 19.230[ms]
u, R = resolver_perthame_gaussiano(
m_1=0.7, m_2=1.1,
eps=0.001, sig_K=0.6, r=1,
M_in=3, mu_in=0, sig_in=1,
u_max=0.6, mu_u=-0.3, sig_u=0.5,
R0_max=0.6, mu_R0=0, sig_R0=1,
solver_u=Solver2U, solver_R=Solver1R,
x_lims=(-4, 4), y_lims=(-4, 4), N=400, M=100, dt=0.01, T=1000
)
plot_total_mass(u, R)
plot_phase_plane(u, R, fig_config={"figsize": (12.5, 5)})
plot_contour_line(u, n_lines=300)
plt.show()
plot_anim(u, R, speed=5.0, t_min=None, t_max=None)
Creating arrays... Building solvers... Starting iterations Iteration n = 999, Total Time: 17.124[seg], Time/Iter: 14.998[ms], avg. Time/Iter: 17.124[ms]
u, R = resolver_perthame_gaussiano(
m_1=0.5, m_2=0.5,
eps=0.001, sig_K=0.05, r=1,
M_in=3, mu_in=0, sig_in=3,
u_max=3, mu_u=-0.3, sig_u=1,
R0_max=5, mu_R0=0, sig_R0=1,
solver_u=Solver2U, solver_R=Solver1R,
x_lims=(-4, 4), y_lims=(-4, 4), N=400, M=100, dt=0.01, T=1000
)
plot_total_mass(u, R)
plot_phase_plane(u, R, fig_config={"figsize": (12.5, 5)})
plot_contour_line(u, n_lines=300)
plt.show()
plot_anim(u, R, speed=5.0, t_min=None, t_max=None, ax_config={"ylim": (0, 4)})
Creating arrays... Building solvers... Starting iterations Iteration n = 999, Total Time: 17.536[seg], Time/Iter: 15.022[ms], avg. Time/Iter: 17.536[ms]